% PSO算法（可以进行调整）
function fitness_history = ae_pso(n_particles, n_variables, max_iter, lb, ub, objective_function)
    
    w_inertia = 0.6;  % 默认惯性权重
    c1 = 0.5;         % 默认个体学习因子
    c2 = 0.5;         % 默认社会学习因子

    N=n_particles;
    D=n_variables;
    MaxFEs=max_iter;
    
    % 粒子位置初始化
    X = lb + (ub - lb) * rand(n_particles, n_variables);  % 粒子位置
    V = rand(n_particles, n_variables);                   % 粒子速度

    % 初始化个体最优解和全局最优解
    Pbest = X;  % 个体最优解
    Pbest_fitness = arrayfun(@(i) objective_function(X(i, :)), 1:n_particles);  % 个体最优适应度
    [gbest_fitness, gbest_idx] = min(Pbest_fitness);  % 全局最优适应度
    Gbest = X(gbest_idx, :);  % 全局最优解
    R=X;
    W=X;
    newE=X;
    f=Pbest_fitness;
    
    % 记录适应度函数的变化
    fitness_history = zeros(max_iter, 1);  % 存储每一代的全局最优适应度值
    Pa=lb+(ub-lb).*rand(1,D);
    Pb=lb+(ub-lb).*rand(1,D);

    % PSO主循环
    for FEs = 1:max_iter
        [~,ind]=sort(f);
        R1=rand(N,D);R2=rand(N,D);S=randi([0,1],N,D);
        r=(ub-lb).*(2*R1.*R2-R2).*S;
        %Computing alpha
        alpha=exp(log(1-FEs/MaxFEs)-(4*(FEs/MaxFEs))^2);
        ar=alpha.*r;
        X=Pbest;
        for i = 1:n_particles
            cab=FEs/MaxFEs;
            %Sampling W and L
            R(i,:)=X(ind(randi([1,length(1:find(i==ind))])),:);
            W(i,:)=X(ind(randi([length(1:find(i==ind)),N])),:);
            if rand<0.5
                A=X(randi(N,D,1),:);
                Pa=(1-cab)*Pa+cab*diag(A)';
                Ov=Pa;
            else
                K=ceil(N*rand);
                I1=[];I1=randperm(N,K);
                w=[];w=f(I1)./sum(f(I1));
                B=[];B=X(I1,:);
                Pb=(1-cab)*Pb+cab*(w*B);
                Ov=Pb;
            end
            I2=round(rand);
            sita=[];sita=I2*rand(1,D)+(1-I2)*rand*2;
            
            %AE Core operator
            newE(i,:)=Ov+ar(i,:)+sita.*(R(i,:)+Pbest(i, :)-Ov-W(i,:));
            Efitness = objective_function(newE(i, :));
            if Efitness < Pbest_fitness(i)
                Pbest(i, :) = newE(i,:);
                Pbest_fitness(i) = Efitness;
            end
                


            % 更新粒子速度和位置
            V(i, :) = w_inertia * V(i, :) + c1 * rand() * (Pbest(i, :) - X(i, :)) + c2 * rand() * (Gbest - X(i, :));
            X(i, :) = X(i, :) + V(i, :);
            
            % 限制粒子位置在搜索空间内
            X(i, :) = max(min(X(i, :), ub), lb);
            
            % 计算新位置的适应度（目标函数值）
            new_fitness = objective_function(X(i, :));

            % 更新个体最优解
            if new_fitness < Pbest_fitness(i)
                Pbest(i, :) = X(i, :);
                Pbest_fitness(i) = new_fitness;
            end
            f=Pbest_fitness;
        end
        
        % 更新全局最优解
        [min_fitness, min_idx] = min(Pbest_fitness);
        if min_fitness < gbest_fitness
            gbest_fitness = min_fitness;
            Gbest = Pbest(min_idx, :);
        end
        
        % 记录每次迭代的全局最优适应度值
        fitness_history(FEs) = gbest_fitness;
    end
end
